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In this paper, low-energy scattering of the ( D*D*) ± meson system is studied within Liischer’s 
finite-size formalism using Nf = 2 twisted mass gauge field configurations. With three different pion 
mass values, the s-wave threshold scattering parameters, namely the scattering length ao and the 
effective range r o, are extracted in ,J F = 1 + channel. Our results indicate that, in this particular 
channel, the interaction between the two vector charmed mesons is weakly repulsive in nature hence 
do not support the possibility of a shallow bound state for the two mesons, at least for the pion 
mass values being studied. This study provides some useful information on the nature of the newly 
discovered resonance-like structure Z c ( 4025) observed in various experiments. 


I. INTRODUCTION 

Since the observation of charged charmonium-like 
structure Z c ( 3900) [T], the BESIII Collaboration stud¬ 
ied the process e + e - -A tt ± (D*D*) ± at a center- 
of-mass energy of 4.26 GeV and reported a new 
charged charmonium-like structure which they named 
as Zf{ 4025) [2], with a mass of 4026.3±2.6±3.7 MeV 
and a width of 24.8±5.6±7.7 MeV. Such charged 
charmounium-like states are quite unique in the sense 
that their valence quark component must contain tetra- 
quark content qiq 2 CC where q± and q 2 being two differ¬ 
ent flavors of light quark. Another feature is that, their 
mass values are rather close to the threshold of two cor¬ 
responding charmed mesons. It is therefore tempting to 
explain these new exotic states as shallow bound states of 
the corresponding mesons. Another explanation is that 
they are simply genuine tetra-quark hadrons or mixture 
of the tetra-quark and the two-meson system. Since it is 
still unclear whether these states are above or below the 
threshold, it is also possible that they are resonances, 
or even simply cusp effects due to interaction between 
different channels. Obviously, a better understanding of 
the internal structures of these states will provide new 
insights into the dynamics of multi-quark systems and 
QCD low-energy behaviors. 

The experimental discovery of the charged 
charmonium-like structures have triggered a lot of 
theoretical studies in recent years, both using phe¬ 
nomenological methods 0HS] and on the lattice [3 M- 
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Since Zf{ 4025) is near the D*D* threshold, a shallow 
bound state, also known as the molecular state, formed 
by D* and D* mesons is a possible explanation. To 
further investigate this possible scenario, the interaction 
between D* and D* mesons at low-energies becomes 
important. The energy considered here is very close to 
the threshold of the D*-D* system. Therefore the inter¬ 
action between the charmed mesons is non-perturbative 
in nature which requires a genuine non-perturbative 
framework such as lattice QCD. 

In this paper, the near-threshold scattering of 
(D*D*) ± system is studied using Nf = 2 twisted mass 
gauge field configurations [9]. The study is carried out 
for three different values of pion mass corresponding to 
m n = 300,420,485 MeV and the size of the lattices is 
32 3 x 64 with a lattice spacing of about 0.067 frn. Accord¬ 
ing to the BESIII results [2], the state 2), (4025) is con¬ 
sistent with the quantum number assignment J p = 1 + 
although other assignments are not completely ruled out. 
Taking into the fact that the state is so close to the 
threshold where presumably s-wave scattering will domi¬ 
nate, we will focus on the J p = 1 + channel only. Exper¬ 
iments also indicates that the state is strongly coupled 
to the D*D* system. Thus, in this exploratory lattice 
study, single-channel scattering of a D* and a D* meson 
is studied using Liischer’s formalism m- The s-wave 
low-energy scattering parameters, namely the scattering 
length ao and the effective range ro, are extracted from 
our simulation. To enhance the energy resolution close 
the threshold, twisted boundary conditions are utilized. 

This paper is organized as follows. In section [TT} we 
briefly recapitulate Liischer’s formalism in general and in 
the particular case of twisted boundary conditions. Sec¬ 
tion |TTT] defines the one-particle and two-particle interpo- 
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lating operators used in this study and the corresponding 
correlation functions. In section |IV[ simulation details 
are provided and the results for the single-meson and 
two-meson systems are analyzed. By applying Liischer’s 
formula, the scattering phases are extracted and when 
fitted to the known low-energy behavior, the threshold 
scattering parameters of the system, i.e. the inverse scat¬ 
tering length Oq 1 and the effective range r 3 are obtained. 
As a crosscheck, both the jackknife and the bootstrap 
method have been used in this study which yield com¬ 
patible results. Implications of our results are discussed 
afterwards. We finally conclude in section [V] with some 
general remarks. 


II. THEORETICAL FRAMEWORK 


Let us first consider a particle with a mass m enclosed 
in a cubic box of size LxLxL, then the ordinary periodic 
boundary condition in the spatial directions reads 

^(x + Lei,t) = \l/(x, t) , (1) 

with the Cartesian unit vector along the i-th axis (* = 
1, 2, 3 for x , y , z direction). The spatial momentum k of 
this particle is quantized according to: 

27T 

k = — li. n € Z 3 . (2) 

L 

Now consider two interacting particles with masses m\ 
and m 2 in this finite box. Taking the center-of-mass 
frame of this system, the two particles thus have oppo¬ 
site three-momentum k and —k. The exact energy E \_2 
of the two-particle system is parameterized as 


-Ea. 2 (k) = \jm\ + k 2 + \Jm\ + k 2 , (3) 


where k 2 is a quantity which also encodes the interaction 
of the two particles in this box. To be specific, k 2 = 
k 2 corresponds to the non-interacting case, while k 2 > 
k 2 and k 2 < k 2 corresponds to repulsive or attractive 
interaction, respectively. Based on Eq. © , it is more 
convenient to define a dimensionless quantity q 2 : 


2 _ k 2 L 2 

q 


(4) 


such that the repulsive and attractive interactions are 
translated into q 2 > n 2 and q 2 < n 2 for some n £ Z 3 , 
respectively. 

In an actual lattice computation, the exact energy 
E\ 2 of the two-particle system hence also the value of 
q 2 is obtained from corresponding correlation functions. 
Liischer’s formula relates the value of q 2 and the elastic 
scattering phase shift S(q) at that particular energy in 
the infinite volume. In the simplest case of s-wave elastic 
scattering, it reads: ilOj 

q cot S 0 (q) = -jj^Z 00 (l;q 2 ) , (5) 


where Zoo(l;q 2 ) is the zeta-function which can be eval¬ 
uated numerically once its argument q 2 is given. Eq. © 
is the main formula to compute the elastic scattering 
phase shift on the lattice. In the case of attractive inter¬ 
action, the lowest two-particle energy level can become 
lower than the threshold. If the interaction is weak, the 
state is loosely bound, i.e. (— q 2 ) being positive but close 
to zero mm- However, a negative q 2 value in a finite 
volume alone does not signifies a bound state. One has 
to investigate the behavior of the negative energy shift in 
the large volume limit. 

With quantization condition on three-momenta, c.f. 
Eq. ©, the typical size of the smallest nonzero momen¬ 
tum is still too large to investigate the hadron-hadron 
near-threshold scattering for practical size of the lattice. 
We thus utilize the so-called twisted boundary conditions 
in our study [13] [14| . Following the notation in Ref. T§] , 
the quark field -^Wx, t ), when transported by an amount 
of L along the spatial direction i (designated by unit vec¬ 
tor e*, i = 1, 2, 3), will acquire an additional phase e zSi : 

-tp 0 {x +Le t ,t) = e l9i i/j 6 (x,t) , (6) 


where 0 = (9i,9 2 ,9 3 ) is the twisted angle (vector) for 
the quark field in three spatial directions. The con¬ 
ventional periodic boundary conditions corresponds to 
6 = (0, 0, 0). Twisted boundary conditions such as those 
in Eq. © can be applied to any flavor of quark fields in 
question. In other words, we are free to choose a twist¬ 
ing angle vector 6j for flavor /, with / = u,d,s,c, - • ■. 
Under twisted boundary conditions, the discretized mo¬ 
mentum in the finite volume is also modified. So, instead 
of Eq. I©, we have, 


P = 


27T 

L 



(7) 


It is more convenient to introduce the new fields '0 / , we 
shall call them the primed fields, via 

^'(x, t) = e~ l6 - x/L tp e {x, t ) . (8) 

It is easy to verify that the primed fields sat¬ 

isfy the usual periodic boundary conditions, c.f. Eq. ©. 
For Wilson-type fermions, we can easily calculate the 
primed quark propagators, which are Wick contractions 
of the primed fields, using a modified set of gauge fields 
(the primed gauge fields), U' x = e * 6 b 0 / i [/ x M with 9 M = 

(o,0)[3[H]. 

Traditional meson interpolating operators are con¬ 
structed using the primed fields as a local bilinears, 
0r(x, t) = ftt) i where / and f denoting fla¬ 
vor indices and T being a Dirac gamma matrix. By 
summing over the spatial coordinate x with appropriate 
three-momentum p, 


C’r(PU) = 5Z^ , / r ^/'( x ’*) e * PX ’ ( 9 ) 

X 


one sees that the above operator in fact corresponds to an 
operator built using the un-primed fields with three mo¬ 
mentum: P+ (Off — 0f)/L. Since it is free to choose any 
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values of Of and Of>, an improved resolution is achieved 
in momentum space. 

Note that we have adopted twisted boundary condi¬ 
tions for the valence quark fields. This is referred to as 
the partial twisting. Strictly speaking, the same twisted 
boundary condition should be applied both to the valence 
and to the sea quark fields which is called full twisting. 
It has been shown recently that, in some cases, partial 
twisting is equivalent to full twisting m- In other cases, 
however, the corrections due to partially twisted bound¬ 
ary conditions are shown to be exponentially suppressed 
if the size of the box is large j!4| . We will assume that 
these corrections are indeed negligible. 0 In the follow¬ 
ing calculations, only the light quark fields (u and d) 
will be twisted while the charm quark fields remain un¬ 
twisted. This choice carefully avoids potential problems 
that might have arisen due to annihilation diagrams in 
this process as suggested in Ref. En¬ 


in. OPERATORS AND CORRELATORS 

As usual, the energies of single-particle and two- 
particle systems are obtained from corresponding corre¬ 
lation functions which are measured in our Monte Carlo 
simulation. Since the newly discovered Z c ( 4025) state 
is observed in both D*D* and the h c tt channel [2], its 
quantum number is likely to be I G (J P ) = 1 + (1 + ). The 
closeness of its mass to the D*D* threshold suggests that 
it might be a candidate for D*-D* bound state. In order 
to investigate the scattering relevant to this scenario on 
the lattice, we need to construct the D*D* two-particle 
interpolating operators with the right quantum number 
mentioned above. In practice, for the one-particle oper¬ 
ators of D * ± and D*°, conventional quark bilinear op¬ 
erators for vector mesons are utilized. The desired two- 
particle operators for system in the I G (J P ) = 1 + (1 + ) 
channel are discussed in the following. Due to the dif¬ 
ference in the symmetries, the cases of twisted boundary 
conditions and non-twisted boundary condition have to 
be treated somewhat differently . 


A. Operators in the non-twisted case 

Let us first consider the non-twisted case. For a single 
vector charmed meson and its anti-particle, we utilize the 
following local interpolating fields in real space: 

[£>*+] : = [c 7 ;d](x,t) , (10) 

[D*~] : Vi{x.,t) = [d 7 ic](x,t) = [T’ l (x,t)] t , (11) 


1 This makes sense since Liischer’s formalism also requires that 
exponentially suppressed corrections are negligible anyway. 


In the above equation, we have also indicated the quark 
flavor content of the operator in front of the definition 
inside the square bracket. So, for example, the oper¬ 
ator in Eq. (10) will create a D* + meson when acting 


on the QCD vacuum. A single-particle state with defi¬ 
nite three-momentum k is defined accordingly via usual 
Fourier transform HZ]: 


—ik-: 


'Pi(kA) = 5 >(x, t ) 

X 

The conjugate of the above operator is: 


( 12 ) 


[7^(k,t)]t = ^[P i (x,t)]t e + Ikx = Pii-Kt) . (13) 


Similarly, for D*° and its anti-particle, we use the follow¬ 
ing operators: 

[D*°] : Qi(x,t) = [cjiu](x,t) , 

[D*°] : Qi(x,t ) = [u7ic](x,t) = [Q i (x,f)] t , 

Qi(k,i) = ^Qi(x,t)e~ lk ' x , (14) 

X 

[Q ? ;(M)]t = ^[Q. i (x,t)]t e + ikx = Qi(-M) . 


For the two-particle operators, in terms of the oper¬ 
ators defined above, we have used the following combi¬ 
nation for a pair of charmed mesons with back-to-back 
momentum, 


P i (k,t)Q j (-k,t)-P i (k,t)Q i (-k,t) , (15) 


with i. j = 1,2,3. On a finite lattice, however, the rota¬ 
tional group SO( 3) is broken down to the cubic group Oh 
and J p = 1 + of the two-particle system is thus reduced 
to T p of the cubic group. To avoid complicated Fierz 
rearrangement terms, we have put the two mesons on 
two neighboring time-slices. Thus, we use the following 
operators to create the state with two charmed mesons, 


OT 1 : 


o Ti 


T x . 


o 


0 3 Tl : 


X) [Pz{R ° k a , t + 1)Q 3 (—i? ok a ,t) 

R£G 

-p 3 (f?ok a ,t + l)Q 2 (-l?ok 0 ,t)] , 

X) [Pi(Rok a ,t+l)Q 3 (-Rok a ,t) 

rgg 

-V 3 {R°k on t+ l)Qi(-Rok a ,t)] , 

E Pi(fiok a ,t+l)Q2(-flok Q ,l) 
R&G 


—ViiR ° k a ,t + l)Qi(—-R o k a , t)] , 

(16) 


where k Q is a chosen three-momentum mode. The index 
a (a = 1, • • • , N) denotes the momentum mode consid¬ 
ered in our calculation. In this particular case, we have 
N = 4. In the above equation, G = Oh designates the 
cubic group and R £ G is an element of the group and we 
have used the notation Rok a to denote the momentum 
obtained from k Q by applying the operation R on k Q . 






4 


Note that in the above constructions, we have not in¬ 
cluded relative orbital angular momentum of the two par¬ 
ticles, i.e. we are only studying the s-wave scattering of 
the two mesons. This is justified for this particular case 
since close to the threshold, the scattering is always dom¬ 
inated by the s-wave contributions. 

B. Operators in the case of twisted boundary 
conditions 

As explained at the end of previous section, we choose 
to apply twisted boundary conditions to the light quarks 
(u and d ) while the charm quark remains un-twisted. 
Single meson operators are the same as in the previous 
subsection except that all the operators are constructed 
using the primed fields. We also set the twisting angle 
for the u and d quark fields to be identical so that their 
lattice propagators are related to each other by a simple 
conjugation in the twisted mass formalism. 

For the two-particle operators, the only difference is 
the discrete version of the rotational symmetry. It has 
been reduced from Oh to one of its subgroups: Ci v , Dih, 
D 2 h, or D^d, depending on the particular choice of 6. The 
other structures (flavor, parity when applicable etc.) of 
the operators remain unchanged. As a consequence, the 
operators V% and Qi, which used to form a basis for the 
Ti irrep of Oh now have to be decomposed into new basis 
of the corresponding subgroups mm- 


Ti Ax ® E 

Civ , 

T\ 1 —y A.2 ® E 

E>ih , 

T\ 1 —y ® B 2 ® B% 

D 2 h , 

T\ 1 —y A .2 ® E 

D 3 d ■ 


The information for these decompositions are summa¬ 
rized in Table [T] As an example, take the first line of 
Eq. @ which corresponds to the case of 0 = ( 0 , 0 , 7 r), 
the original operator triplet (V\ ,V2,V 3 ) should be decom¬ 
posed into a singlet ( V 3 ) and a doublet (V i,V 2 ) which 
forms the basis for A\ and E irreps, respectively. Similar 
relations also hold for the Qi s. 

The construction of the two-particle operators in the 
case of twisted boundary conditions is somewhat com¬ 
plex. Let us start from a general problem in group theory. 
Suppose that ei form the basis of a 3-dimensional irreps 
Ti while e[ form the basis of another 3-dinrensional irreps 
T\. With the help of group theory, the direct product of 
these two 3-dinrensional irreps can form a 9-dimensional 
reducible representation of basis a g> e'(*,j = 1,2,3). 
Depending on the particular choice of 0, this new 9- 
dimensional reducible representation will be decomposed 
into irreps of the corresponding subgroup, with the linear 


2 The reason that V 3 is special as opposed to P\ and P 2 is because 
the twisted boundary condition with 6 = (0,0, n) is applied in 
the 3-direction which breaks the symmetry. 


combinations of e* 0 e'(*, j = 1,2,3) giving the basis of 
these irreps. 

To find the linear combination of basis ej (g> e(-(i, j = 
1,2,3) for definite irrep we are interested in, one could 
use different approaches. In our study, group character 
technique is used to determine the specific basis for a 
certain irrep. 

As an application of this technique described above, 
taking the case of Q = (7r, 7 r, 0) as an example, we give the 
corresponding operators as listed in the following equa¬ 
tions. 

Bi : ei g> e' 2 + e 2 (g> e[ , 

B 2 : ei g> e 3 + e 3 (g> e[ , (18) 

B 3 : e 2 ® e 3 + e 3 <g> e' 2 , 

where d = ^ 75(^1 + V 2 ), e 2 = 75(^2 ^ Vi), e 3 = V 3 . 

Similar relations also hold between e! i and Qi. Then we 
have two-particle operators for irrep B 1 as shown below: 

0 Bl : [P 2 (R o k a , t + 1)Q 2 (—R o k a , t) 

R&G (19) 

— Vi(Ro k a , t + l)Qi(— 1? o k a , f)] . 

where G = D 2 h, the group corresponding to 0 = ( 7 r, 7 r, 0 ). 

C. Correlation functions 

For vector charmed meson D* and D* , the correspond¬ 
ing correlation functions are defined as: 

C v (Kt) = (v}(Kt)Vi(K 0 )), 

C c (k,t) = (Ql(k,t)Q l (k,0)) , 

where k represents the three-momentum of the relevant 
particle. It is straightforward to obtain the single particle 
energy E( k) for various lattice momentum k. For the sin¬ 
gle particle, the dispersion relation can then be checked 
with various E{ k). In particular, this can be checked in 
both twisted boundary conditions and conventional pe¬ 
riodic boundary conditions. With judicious choices of 6, 
one could check the single-particle dispersion relation to 
a much better accuracy which will be shown in the next 
section. 

Two-particle correlation functions are somewhat more 
involved. Generally speaking, a correlation matrix 
Cl p {t) is constructed: 

Clp{t) = (O r J(t)O B (0)) . (21) 

where O £ represents the two-particle operator defined 
in the previous section and T denotes a definite irrep 
while a enumerates different operators in that irrep. To 
be specific, for the non-twisted case 0 = ( 0 , 0 , 0 ), the 
number of k a is 4 in Ti channel while for all other cases, 
the number of k„ is 2 . As a reference, these information 
are also collected in Table |U 
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TABLE I. Information about the two-particle operators used 
in this calculation together with the corresponding symme¬ 
tries. 


e 

0 

(0,0,?) 

(0,0,7r) 

(7T, 7T, 0) 

(7T, 7T, 7t) 

Symmetry 

o h 

C^V 

D^h 

D 2 h 

D^d 

irreps 

T, 

Ai, E 

E 

B\, B 2 , B 3 

E 

Number of k,, 

4 

2 , 2 

2 

2,2,2 

2 


IV. SIMULATION DETAILS AND RESULTS 

In this paper, the Osterwalder-Seiler action m is used 
for the valence charm quark. The gauge field ensemble 
comes from Nf = 2 twisted mass gauge field configura¬ 
tions generated by the European Twisted Mass Collabo¬ 
ration (ETMC) jSJ. The gauge coupling is j3 = 4.05 which 
corresponds to a lattice spacing of about 0.067 fm and 
we have used three different pion mass values, namely 
300 MeV, 420 MeV and 485 MeV. Details of the rele¬ 
vant parameters are summarized in the Table [TTJ The up 
and down bare quark mass values, characterized by the 
bare quark parameter p in Table [TIJ are fixed to that of 
the sea-quark. For the charm quark, the mass parameter 
a/j, c is fixed so that the value of \m Vc . +1 mj/^ calculated 
on the lattice reproduces the corresponding experimental 
value. 


TABLE II. Simulation parameters in this study. 



-^Vconf 

rriTr [MeV] 

m n L 

Id 

X 

T 

a[fm] 

P 

0.003 

200 

300 

3.3 

3 2 3 

X 

64 

0.067 

4.05 

0.006 

200 

420 

4.6 

32 3 

X 

64 

0.067 

4.05 

0.008 

200 

485 

5.3 

32 3 

X 

64 

0.067 

4.05 


A. One-particle spectrum and dispersion relation 

One-particle correlation functions as defined in 
Eq. ( [20| with definite three-momentum k are calculated 
in our simulation from which the one-particle spectrum 
E(k) is obtained. We have checked the single particle 
dispersion relations for D* and D* mesons, with both 
periodic boundary conditions and twisted boundary con¬ 
ditions. For the twisted boundary conditions, equivalent 
small momentum points offer us a more stringent test for 
the dispersion relation, both the continuum one and its 
lattice counterpart, at low-momenta close to zero. One 
example of these is illustrated in Fig.[l]at p, = 0.003. The 
quantity E( k ) 2 or its lattice counterpart 4sinh 2 (S/2) 
is shown versus p 2 or p 2 = 4 JT sin 2 (pj/2) in the bot¬ 
tom/top panel, respectively. The straight lines are linear 
fits with Z being the fitted slope of the lines. 


dispersion (p=0.003) 



dispersion (p=0.003) 



FIG. 1. Dispersion relation for the D* meson at p = 0.003 
with lattice case (upper panel) and continuous case (lower 
panel). The points with error bars are lattice data while the 
straight lines are the corresponding linear fits with Z denoting 
the slope of the line. 

B. Extraction of two-particle energy levels 

In this paper, the usual Liischer-Wolff method [T9] is 
adopted to extract the two-particle energy eigenvalues. 
For this purpose, a new matrix fl(f, to) is constructed as: 

n(t,t 0 ) = C(to)-*C(t)C(to)~i , ( 22 ) 

where to is a reference time-slice. Normally to is picked 
such that the signal is good and stable. The energy eigen¬ 
values for the two-particle system are then obtained by 
diagonalizing the matrix fl(f,f 0 )- The eigenvalues of the 
matrix, A a (t,to), have the usual exponential decay be¬ 
havior as described by A a ~ e -E a (t-t 0 ) anc [ therefore the 
exact energy E a can be extracted from the effective mass 
plateau of the eigenvalue A a . 

The real signal for the eigenvalue in our simulation 
turns out to be somewhat noisy. To enhance the signal, 
the following ratio was attempted: 

7? M- ^a(t, t 0 ) 

a[,0> C7^(t-to,0)CG(t-t 0 ,0) (23) 

OC e - AE °‘'( t ~ t o) ; 

where C v (t — to,0) and C®(t — to,0) are one-particle 
correlation functions with zero momentum for the corre¬ 
sponding mesons defined in Eq. ( |~i~2j ) and Eq. ( |T4| ). There¬ 
fore, A E a is the difference of the two-particle energy 
measured from the threshold of the two mesons: 


A E a = E a — raj). — m^. ■ 


(24) 
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The energy difference A E a can be extracted from the 
plateau behavior of the effective mass function A E a ^ e s(t) 
constructed from the ratio TZ a (t. t u ) as usual: 


A E_D 4h _E (fj. = 0.008) 


A£ , Qjeff (t) = In ^ 


K Q (t,to) 

lZ a (t + 1, to) 


(25) 


With the energy effective energy difference A E a ,eff(i) 
for each time slice t, we estimate the error for each 
AB 0 ,eff(t) using the jackknife method. Then, from the 
effective energy difference A E ate g(t) and its correspond¬ 
ing errors, one searches a plateau in t that extends several 
consecutive time-slices and minimizes the y 2 per degree 
of freedom. From this procedure, a fitted value of A E a 
together with its error is obtained. As an illustration, 
in Fig. [2j we have shown the fitted values of A E a using 
the jackknife method at fj, = 0.008 in the E channel for 
9 = (0,0,7r), and the B 2 channel for 0 = ( 77 , 77 , 0). As a 
cross check, bootstrap method is also tried to calculate 
the standard error of A E a ^g(t) on each time slice. To 
make the comparison, the fitted values of A E a are also 
illustrated in Fig. [3] for the same cases as in Fig. [5] 




AE_D^j l _E (/i — 0.008) 




FIG. 2. Effective mass plots for the energy shift A E a us¬ 
ing the jackknife method at /j, = 0.008 in the E channel for 
6 = (0,0,7 r) (top) and B 2 channel for 9 = ( 77 , 77 , 0) (bottom). 
The red crosses and the blue open circles correspond to two 
different energy levels obtained from Eq. ( |2 1 1 ) using N = 2 
different two-particle operators as discussed in the text. The 
horizontal bands indicate the fitted values for AE n and the 
corresponding fitting ranges. 

It is seen graphically from Fig. [2] and Fig. [3j the fitted 
values of A E a from jackknife method is consistent with 
those from bootstrap method within the statistical un¬ 
certainties. The only difference is that bootstrap method 
seems to give a somewhat smaller error of A E a ( t ) on each 


FIG. 3. Same as Fig. [2] but using the bootstrap method. 

time slice. To be on the safe side, in this paper we regard 
the results from jackknife method as our final results for 
the energy levels. 

Effective mass plots for other cases are similar to those 
shown in Fig. [2] and Fig. [3j With the energy difference 
AE a extracted from the simulation data, one utilizes the 
definition: 

\J m 2 D . + k 2 + yj m 2 B , + k 2 = AE a + rriD* Amp* • (26) 

to solve for k 2 = ( 2n/L) 2 q 2 which is then plugged into 
Liischer’s formula to obtain the information about the 
scattering phase shift. 

The final results for A E a in each irrep, together with 
the corresponding ranges from which the AE a ’s are ex¬ 
tracted, are summarized in Table m We only keep the 
two lowest energy levels for the non-twisted case and the 
lowest for the twisted cases, since those higher energy lev¬ 
els are not going to be utilized to extract the scattering 
parameters in the following analysis anyway. |^] As a re¬ 
sult, altogether 9 energy levels are kept for the scattering 
analysis in the following. 

C. Extraction of scattering information 

The energy considered in this study is very close to 
the threshold of the D*-D* system, therefore one has the 


! The additional operators in each particular channel helps to sta¬ 
bilize the lowest energy values in each irrep although the actual 
values of these higher states are not utilized. 
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e 

Irrep 

AS[t min , £ max ](/u = 0.003) 

A S[t min ,t max ](/x = 0.006) 

^max 

](/z = 0.008) 



pmodeO 

pmodel 

pmodeO 

pmodel 

pmodeO 

pmodel 

0 

Ti 

0.001(2) [8,13] 

0.068(3) [6,11] 

0.005(2) [9,14] 

0.059(4) [8,13] 

0.003(1) [8,13] 

0.046(3) [8,13] 

fc\ n 77 \ 

Ai 

0.001(3) [8,15] 


0.012(2) [6,11] 


0.007(2) [8,13] 


u, 2 j 

E 

0.008(2) [6,11] 


0.013(1) [6,11] 


0.008(2) [9,14] 


(0,0, tt) 

E 

0.004(5) [10,15] 


0.019(2) [8,13] 


0.017(2) [9,14] 



Si 

0.027(3) [8,15] 


0.035(3) [9,15] 


0.027(1) [7,12] 


(7T, 7T, 0) 

s 2 

0.032(3) [8,13] 


0.030(3) [9,15] 


0.025(2) [8,14] 



S 3 

0.032(2) [7,13] 


0.038(2) [7,12] 


0.028(2) [9,14] 


(tt, TT, 7t) 

E 

0.040(2) [5,10] 


0.047(1) [5,10] 


0.038(3) [7,12] 



TABLE III. Results for the energy shifts A E obtained in our calculations for various cases. The time interval [i m i n ,imax] from 
which we extract the values of A E are also listed. These ranges are relevant for the estimation of the error for the zeta functions 
as described in the text. 


following effective range expansion: 


1 


k 2l+1 cot Si(k) = 1 + -rik 2 + 


(27) 


where ai is the so-called scattering length for partial wave 
l and ri is the corresponding effective force range while 
• • • represents terms that are higher order in k 2 . It is 
more convenient to use a dimensio nless form in our anal¬ 
ysis. With q 2 = k 2 L 2 /( 27r) 2 , Eq. (27) can be rewritten 
in terms of q 2 : 


This yields the low-energy scattering parameters 
for both s-wave and p- wave; 

2. A correlated fit using only the parity-conserving 
points. This yields only the s-wave low-energy scat¬ 
tering parameters. 

3. A correlated fit using all data points, neglecting 
the parity-mixing effects of the two data points for 
9 = (0, 0,7t/ 2). This also only yields the s-wave 
scattering parameters. 


q 2l+1 cot Si(q 2 ) = Bi+ * Riq 2 -f-, (28) 

with Bi = [L/{2 'k)] 21+1 a^ 1 and Ri = [L/(2'ir)] 2l ~ 1 ri. In 
the following, we will call parameters Bi and Ri the low- 
energy scattering parameters in partial wave l and our 
task is to extract these parameters from our simulation 
data. Since we have a definite lattice size and lattice 
spacing, it turns out that q 2 ~ 1 corresponds to k 2 ~ 
[580MeV] 2 in physical units. 

It is also well-known that, close to the threshold, scat¬ 
tering is dominated by phase shifts coming from lower 
partial waves as long as they are non-vanishing. There¬ 
fore all partial waves with l > 2 will be ignored in the 
Liischer formula for this study. As mentioned in pre¬ 
vious section, the irreps studied in this paper all pre¬ 
serve parity except for the case of 0 = (0,0,7r/2) which 
breaks parity. Using the terminology in Ref. [3], this is 
the only parity-mixing scenario while all other points be¬ 
long to the parity-conserving scenario. Thus to extract 
these low-energy scattering parameters from the lattice 
data, we have altogether 9 points for different q 2 values: 
2 points in the parity-mixing case with 9 = (0, 0, 7t/ 2) 
and 7 points in the parity-conserving case. These are all 
tabulated in Table Q3H 

As all contributions from l > 2 partial waves have been 
neglected, the parity-conserving data (7 points) will de¬ 
pend only on the s-wave parameters Bq and i?o while the 
parity-mixing data (2 points) will depend on both the s- 
wave parameters and the p -wave parameters B\ and Ri . 
There are 3 different strategies to follow here: 

1. A combined correlated fit using all 9 data points. 


We will first describe the fitting process following strat¬ 
egy 1 listed above. The other strategies follow similarly 
and the results will also be listed for comparisons. 

To be specific, in the parity-conserving case, we define 

yo(q 2 ) = m. 00 (q 2 ) . (29) 

According to Liischer’s formula Eq. ([5]), this should be 
equal to 

q cot 5 0 (q 2 ) = m 00 (q 2 ) = -^Z 00 (l; q 2 ) , (30) 

for the non-twisted case while for the twisted case, 
one simply replace the corresponding zeta function by 
Zq Q ( 1; q 2 ) [153. In the parity-mixing case, however, things 
are more complicated. Apart from the s-wave phase shift 
5o(q 2 ), Liischer formula will also involve 8\(q 2 ). Accord¬ 
ingly, we define 

yi(q 2 ) = [moi {q 2 )] 2 , (31) 

and, according to Liischer’s formula, it is equal to 

yi(q 2 ) = [qcots 0 {q 2 ) - m 00 ] [<7 cot <5i (g 2 ) - mu) (32) 

where the functions moo, moi and mu are related to the 
corresponding zeta-functions, see e.g. Ref. m- 

For definiteness, we label the data points as follows: 
the parity-conserving data points are labelled from 1 to 
Nq = 7 while the parity-mixing points are labelled from 
-/Vo + 1 = 8 to N 0 + Ni = 7+2 = 9. For later convenience, 
we also introduce an index function as follows, 

0 for 1 < / < Nq 
1 for 7V 0 + 1 < / < N 0 + IV] 


( 33 ) 





















y = 0.003 


In other words, ind(I) = 0 for the first Nq parity- 
conserving data points while ind(I) = 1 for the next 
Ni parity-mixing data points. So our previous defini¬ 
tions of yo{q 2 ) and yi(q 2 ) may be written collectively as 
Vind(i){qj ) with I = 1, 2, • • • , (N 0 + iVi). We can then 
construct the y 2 function as usual 


No + Ni 

x 2 = £ [Find(i)(qh a) - yind(i)(qj)\ cj} 

i , j -i 

\Find(J){,qj i — Vind(J){qj)^ ■ 


(34) 


where for ind(I) =0,1 the corresponding functions are 
(using the symbol a to collectively denote all the relevant 
fitting parameters I?o, -Ren Ri and Ri): 


F 0 (q 2 ;a) = B 0 + ^q 2 


(35) 


Rq 2 l rn , Rl 


Fi(q ;a) = [B 0 + —q -m 00 ][B 1 + —q -mn].(36) 

For the estimation of the covariance matrix Cjj and also 
the errors for the zeta-functions that appear in the above 
formulas, we closely follow the steps outlined in Ref. iH]. 
The reader is referred to that reference for further details. 


Basically, minimizing the target y 2 function in Eq. (341, 


one could obtain all the parameters, namely Bq, Rq, Ri 
and R\, in a single step with all of our data. In the 
course of inverting the covariance matrix C in Eq. |34[ 
the eigenvalues of the covariance matrix for each irrep 
have been calculated with both QR decomposition and 
singular value decomposition. The corresponding results 
show that the matrices are nonsingular. 

To get a feeling of these fits, we plot the quantity 
< 7 cot S 0 (q’ 2 ) vs. q 2 in Fig. [^obtained from strategy 1. The 
values of <7 cot <5o(<? 2 ) for the data points are obtained via 
the relation 


<7cot<5 0 (<7 2 ) = m 00 


m\ 


01 


q 3 cot Si(q 2 ) — mu 


(37) 


where the quantity q 3 cot<5i(g 2 ) on the r.h.s of the above 
equation is replaced by B\ + \R\q 2 with the fitted values 
for B\ and R\. This figure illustrates the situation for 
all three pion masses in our simulation. From top to bot¬ 
tom, each panel corresponds to y = 0.003, n = 0.006 and 
/r = 0.008, respectively. All data points obtained from 
our simulation are plotted in these figures. The blue open 
circles are the data points in the parity-conserving cases 
while the two red crosses in each panel are the data for 
the parity-mixing case. The straight lines in the figure 
illustrates the fitting function F 0 (q 2 - : a) = B 0 + (R 0 /2)q 2 
and the shaded bands indicate the corresponding uncer¬ 
tainties. As is seen from the figure, we do get a reason¬ 
able fit for all three pion mass values. Finally, the fitted 
values for the scattering parameters are summarized in 
Table IV for three values of m 2 in our simulation. 


To check the validity of the effective range expansion, 
we may also compare the s-wave phase shift 5o(q 2 ) itself 



FIG. 4. Results for the correlated fits from strategy 1 as 
described in the text. Each panel, from top to bottom, corre¬ 
sponds to y, = 0.003,0.006 and 0.008, respectively. The quan¬ 
tity qcotSo{q 2 ) is plotted versus q 2 for all our data points, 
both parity-conserving case (blue open circles) and parity¬ 
mixing case (red crosses). The straight lines indicates the 
fitted result for qcot5o(q 2 ) = Bo + (Ro/2)q 2 and the shaded 
bands indicates the corresponding uncertainties. 



Bo 

Ro 

B\ 

Ri 

X 2 1dof 

0.003 

-0.47(35) 

-0.051(243) 

0.32(17) 

-9.15(2.46) 

1.66/5 

0.006 

-0.46(23) 

-1.46(1.38) 

-0.14(05) 

-0.59(32) 

9.92/5 

0.008 

-0.83(17) 

3.18(2.12) 

0.85(31) 

-14.42(5.83) 

3.52/5 


TABLE IV. Fit results with strategy 1. 



B 0 

Ro 

X 2 /dof 

0.003 

-0.42(38) 

-0.13(43) 

1.62/5 

0.006 

-0.47(19) 

-1.45(1.16) 

9.92/5 

0.008 

-0.84(13) 

3.18(2.29) 

3.52/5 


TABLE V. Fit results with strategy 2. 


as a function q 2 . 


The situation is shown in FIG. [ 5 ] for 
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Bo 

Ro 

x'/dof 

0.003 

-0.57(27) 

0.15(59) 

2.44/7 

0.006 

-0.33(21) 

-1.72(1.25) 

10.94/7 

0.008 

-0.61(3) 

2.94(2.58) 

5.04/7 


TABLE VI. Fit results with strategy 3. 


p = 0.006 



—0.72(54)fm, -0.74(37)fm, -0.41(8)fm for p = 0.003, 
0.006, 0.008, respectively. The values for ro can be also 
obtained accordingly. These numbers are summarized in 
Table Em 

From strategy 2 and strategy 3, the correlated fitting 
are also conducted to obtain the scattering length ao and 
effective force range 7'o- The specific results are summa¬ 
rized in Table |VTIT| and Table |IX| respectively. 



p = 0.003 

p = 0.006 

p — 0.008 

a o [fm] 

-0.72(54) 

-0.74(37) 

-0.41(8) 

r 0 [fm] 

-0.018(83) 

-0.50(43) 

1.08(73) 


TABLE VII. The values for a o and ro in physical units ob¬ 
tained from the numbers for the correlated fit in Table m 
from strategy 1. 


FIG. 5. The same as Fig. [4] but the comparison is done for the 
s-wave phase shift 5o(q 2 ) itself. This is the case of /r = 0.006. 


p = 0.006. 

Similarly, one could follow strategy 2 listed above 
and obtain the s-wave parameters using only the parity- 
conserving data points, i.e. the first 7 data points. Or, 
alternatively following strategy 3 and obtain the s-wave 
scattering parameters using all the data points by ne¬ 
glecting the mixing between the s-wave and p-wave. Nu¬ 
merically this amounts to setting the matrix elements 
moi = 0 compared with the diagonal ones. The results 
obtained from these strategies can be compared with 
what we get from strategy 1. It turns out that, the mix¬ 
ing of the s- and p -wave indeed has little impact on the 
final results for the s-wave scattering parameters. The 
corresponding fitted results are summarized in Table |V| 
and Table [VT| respectively. 

As is seen from these tables, as far as the s-wave scat¬ 
tering parameters are concerned, it seems that the parity- 
conserving data dominate the final fitting results. This is 
illustrated by consistent values for B 0 and R 0 in Table [TV] 
and Table [V] Finally, we regard our correlated fits from 
strategy 1 with all of our data as being more reliable and 
they are taken as our final results for this paper. 


D. Physical values for the scattering parameters 

The relation between the fitted values of Bi , Ri for 
l = 0,1 and scattering parameters can be expressed as 
follows: 


ai = 



n = Ri 



(38) 


Taking the numbers of correlated fitting in Table |IV[ 
for the s-wave, we can obtain the scattering length oq: 



p = 0.003 

p = 0.006 

p = 0.008 

a 0 [fm] 

-0.80(71) 

-0.73(30) 

-0.41(6) 

r 0 [fm] 

-0.043(146) 

-0.49(46) 

1.09(78) 


TABLE VIII. The values for a o and ro in physical units ob¬ 
tained from the numbers for the correlated fit in Table [V] 
from strategy 2. 



p = 0.003 

p = 0.006 

p = 0.008 

a 0 [fm] 

-0.60(28) 

-1.03(65) 

-0.56(3) 

ro [fm] 

0.05(20) 

-0.59(43) 

1.00(88) 


TABLE IX. The values for ao and r o in physical units obtained 
from the numbers for the correlated fit in Table m from 
strategy 3. 


E. Scattering parameters using the bootstrap 
method 

The errors used in the analysis discussed so far are es¬ 
timated using the jackknife method. To crosscheck these 
results, bootstrap method is also utilized to analyze di¬ 
rectly the final scattering length ao and effective range 
ro- The specific procedure is as follows. 

1. Select randomly 200 configurations from the given 
configurations in Table [TT] for each parameter p. Do 
the selection N random times and label each sample 
by an integer i, i = 1, 2, • • • , N random . 

2. For each randomly selected sample i, repeat the 
analysis process described so far in section [TV] This 
yields one set of scattering parameters, say l/og 
and T-g. 

3. Analyze the distribution of these values. Taking r 0 
as an example, find the values p and q so that these 
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bracket the central 68% of the r l 0 values: 


N(r l 0 < p) 


N, 


= 0.16 


N(r' 0 > q) 


random 


N, 


= 0.16 


random 


(39) 


where N(tq < p) denotes the number of satisfy¬ 
ing r l 0 < p. 


4. the bootstrap estimate of the asymmetric errors for 
the quantity can be given as: 


ro = (r 0 ) 


+(g-( r o)) 

-((ro)-p) 


(40) 


with (ro) denoting the weighted mean of {tq}. 

In this work, we take N ran dom = 60 at three different /i 
(/x = 0.003,0.006,0.008) to estimate the bootstrap error 
of scattering parameters with strategy 2 described in sec¬ 
tion IV. As an illustration, the distribution of 1 /a l 0 and Tq 
for case p = 0.003 are shown in FIG. [6] Meanwhile, the 
asymmetric error of ao and r o are estimated with these 
60 samples. These final specific values are summarized 
in Table IXl 

Alternatively, we have also estimated the bootstrap er¬ 
ror of scattering parameters with strategy 3. The only 
difference is that the data points corresponding to Ai 
and E irreps (6 = (0, 0, 7t/2)) are not left out during the 
process, that is neglecting the parity-mixing of effects of 
the two data points. The final results are summarized in 
Table El 


20 

18 

16 

14 

12 

feao 

8 

6 

4 

2 


-2 -1.5 -1 -0.5 0 

l/aol/m' 1 ] 


fi = 0.003 



/i = 0.003 



FIG. 6. Distribution of {^} and {tq} from strategy 3 for 
a o 

bootstrap method at / j , = 0.003. 

As observed from Table [iXj Table [VniJ Table El and 
Table [Xj jackknife method and bootstrap method yield 
compatible results. 



p = 0.003 

p = 0.006 

p = 0.008 

a 0 [fm] 

—0 74+ 0 - 21 
U - '^-0.15 

—0.68i°;)i 

—0.48±g;H 

ro [fin] 

-0.0048t°'i® 

—0.038±g;i| 

PI 7*3+0-58 
u - * °-0.62 


TABLE X. The values for oo and ro from strategy 2 with 
parity-conserving data points with bootstrap method at p = 
0.003,0.006,0.008. 



p = 0.003 

p = 0.006 

p = 0.008 

a 0 [fm] 

-0.76io.2i 

-O. 86 io.l 2 

-0.59io.li5 

ro [fm] 

— 0.0022ig')g 

-0.14i°;l* 

0.64l°; 5 5 ? 


TABLE XI. The values for a o and ro from strategy 3 with all 
data points with bootstrap method at p = 0.003, 0.006,0.008. 


F. Implication of our results 

As is said, we take the fitted result from strategy 1 as 
our final results, i.e. those in Table |IV| and Table |VII| 
Based on our results, the values of ao do not seem to 
follow a regular chiral extrapolation pattern, at least not 
within the range that we have studied. We therefore kept 
the individual values for ao and ro for each case. This 
irregularity might be caused by the smallness of the value 
m^L ~ 3.3 for p = 0.003. To circumvent this, one has to 
study a larger lattice. 

The negative values of the parameter Bq (hence the 
scattering length ao) indicates that the two constituent 
mesons for the ( D*D*) ± system have weak repulsive in¬ 
teractions at low energies. Therefore, our result does not 
support the bound state scenario for these two mesons. 

Another check for the possible bound state would be to 
look for those negative q 2 values we obtained which cor¬ 
responds to the negative values of 5E listed in Table m 
However, the q 2 for different channel in our study are 
all positive which contradicts the possibility of a bound 
state. Since the cases we are studying is still far from 
the physical pion mass case, we still cannot rule out the 
possibility the appearance of a bound state once the pion 
mass is lowered (and the lattice size L is also increased ac¬ 
cordingly to control the finite volume corrections). Such 
scenarios do occur in lattice studies of two nucleons. 


V. CONCLUSIONS 

In this paper, the low-energy scattering of D* and D* 
is studied with Nf = 2 twisted mass fermion configura¬ 
tions. In our calculation, three different pion mass values 
(m n = 300,420,485 MeV) are utilized to investigate the 
pion mass dependence, and the corresponding lattice size 
is 32 3 x 64 with a lattice spacing a ~ 0.067 fm. We have 
used twisted boundary conditions to enhance the momen¬ 
tum resolution close to the threshold. Using Liischer’s 
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finite-size technique, the s-wave scattering in the chan¬ 
nel J p = 1 + is studied and the scattering parameters 
are obtained by correlated fitting procedure. As a cross¬ 
check, two different statistical error estimating methods, 
jackknife and bootstrap method are utilized which yield 
compatible results for the s-wave scattering parameters. 
The results from a correlated fit with all of the data with 
errors estimated using the jackknife method is regarded 
as the final result for this paper. 

Our results indicate that, for all three pion mass val¬ 
ues that we simulated, the scattering lengths are negative 
which indicates a weak repulsive interaction between the 
the two mesons {D* and D* or its conjugated systems un¬ 
der C-parity or G-parity). Thus a bound state of the two 
mesons in J p = 1 + channel is not supported based on 
our current lattice results. However, as we pointed out al¬ 
ready, we cannot rule out the possibility of a bound state 
for the two vector charmed mesons when the pion mass 
is lowered and the volume is increased accordingly. This 
requires further systematic lattice studies. Furthermore, 
it is also possible that more complete set of interpolation 
operators and a coupled channel study is required. In 
summary, this lattice study has shed some light on the 


nature of (4025) however it remains to be clarified by 
future more systematic studies. 
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